%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Reproducing Table 4 of the paper:
%Reis, Ricardo (2006) "Inattentive Producers," Review of Economic Studies.
%Written by: Ricardo Reis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clc

%  STEP 0: Setting parameters
p=2;					%Lag length in VAR
final=167;              %last date to include in estimation sample
                        %167 is 1995Q4, 183 is 1999Q4
if final==183
    display('Estimate from 2000q1 onwards')
elseif final==167
    display('Estimate from 1996q1 onwards')
end
lambda=0.25;
alpha=0.11;
upper=100;

%  STEP 1: LOADING DATA

%1st column has index, 2nd has dates, 3rd has m, 4th has p, 5th has a, 6th has g_m, 7th has pi, 8th has g_a
load InPro_data.txt
index=InPro_data(:,1);
y=InPro_data(:,3)-InPro_data(:,4);
a=InPro_data(:,5);
gap=y-a;
temp=length(gap);
dgap=gap(2:temp,1)-gap(1:temp-1,1);
pi=InPro_data(:,7);
dates=InPro_data(:,2);
% Data for VAR which comes as a Tx4 vector, with:
% - inflation in the first column
% - HP-filtered real output in the second
% - dates in the third column
VARdata=[pi(2:temp) dgap dates(2:temp)];
FULLdata=VARdata;
N=length(FULLdata);
VARdata=VARdata(1:final,:);         %Reducing the sample to only go until the end of estimation period
T=length(VARdata);		            %Setting number of observations
TT=T-p;					            %Estimation sample length
clear InPro_data10 y a pi dates temp

%  STEP 2: RUNNING THE VAR

%Follow Hamilton, pages 292-293 to run the regression
%The y_t is just row(t) in the data matrix
%Stack all the y_t vertically as columns in matrix Y
Y=zeros(2,TT);				        %initialize matix, to ensure right dimension
for i=1:TT
    Y(1:2,i)=VARdata(T+1-i,1:2)';	%first column has most recent observation, so I am going from T to p
end
clear i
%The x_t are 1+2*p vectors, described in eq. 11.1.5 in Hamilton page 292
%These are stacked horizontally to form matrix X which is  1+2*p x TT
X=zeros(2*p,TT);			%initialize matrix
X=[ones(1,TT);X];			%add column of 1's to take account of the constant
for j=1:p					%loop over the different lags
    temp1=2+(j-1)*2;		%index on block will get
    temp2=temp1+1;
    for i=1:TT
        X(temp1:temp2,i)=VARdata(T+1-i-j,1:2)';	%loop over the different dates
    end
end
clear temp1 temp2 j i
%Calculate PI, the matrix of estimated coefficients, using eq. 11.1.11 in
%Hamilton page 293.
sumy=zeros(2,1+2*p);			%sumx will have the sum of x*x' over T
sumx=zeros(1+2*p,1+2*p);	    %sumy will have the sum of y*x' over T
for i=1:TT						%iterate over the estimation sample
    xt=X(:,i);					%xt is the t'th obs. on the x's
    xx=xt*(xt');				%this is the desired product
    sumx=sumx+xx;				%and this adds it up to form the sum
    yt=Y(:,i);					%Get the t'th obsv on the y's then
    yx=yt*(xt');				%do the multiplication
    sumy=sumy+yx;				%and add it up as well
end
PI=sumy*inv(sumx);			    %then these are the estimated coefficients
PI=PI';                         %Note that it is the transpose of Hamilton eq 11.1.6
clear xi xx sumx yi yx sumy xt yt i Y X 
%PI
%display('This is the estimated VAR. Check this should be a (2+2*p) x 2 matrix.' )
%display('Per equation there were 2+2*p regressors, and there are two equations.')
%display('See equation 11.1.6 in Hamilton for what each coefficient refers to.')

%STEP 3: GENERATING FORECASTS

%a) Rewriting the VAR(p) as a VAR(1)
%To do so I will use the trick of re-writing a VAR(p) as a VAR(1)
%which involves building an F matrix described in Hamilton, page 259, eq. 10.1.10
F=zeros(2*p,2*p);					%This matrix is of size 2*p x 2*p
F(1:2,1:p*2)=PI(2:2*p+1,1:2)';	    %Its first block row has the PI without the constant
temp=eye(2);						%this creates the identity matrix of size 3x3
for i=1:p-1							%iteration to put the identity matrices in place
    temp1=1+i*2;					%starting location
    temp2=temp1+1;					%finishing location
    F(temp1:temp2,temp1-2:temp2-2)=temp;
end
clear temp temp1 temp2 i
%F
%display('Check if this F matrix has the form of Hamilton, page 259, 10.1.10')
%b) Building matrix with forecasts
%The matrix Forecasts has the matrix with the forecasts.
%Column 1 has the forecasts of y_t as of t-j+1, where j stands for the row.
%Column 2 has the forecasts of y_t-1 as of t-j, where j is the row.
%The T element is just the constant in the regression, since for y_t at date -1 is
%just like the unconditional expectation since have no observations to condition on.
temp=eye(2);
for i=1:p
    temp=temp-PI(2*i:2*i+1,1:2)';            %calculate the unconditional forecast, mu, following Hamilton, p. 258 
end
mu=inv(temp)*PI(1,1:2)';
clear temp i
Forecasts=[];
for t=1:N				%iterate over the y_t the variable trying to forecast
    j=N-t+1;			    %set it up so that t=1 means the last observation.
    piForecast_t_temp=[];
    dgForecast_t_temp=[];
    for s=1:N							%indexing over the period of forecast: y_t_t-s
        pi_for=0;                       %initializing the final forecast
        yt_s = mu;                      %adding constant
        F_forecasting=F^(s-1);          %this is what allows for forecasting k periods into the future
            for i=1:p					%add the lagged observations to the forecast
                if j-s+2-i>0		    %but only if they exist in the sample
                    yt_s=yt_s+F_forecasting(1:2,2*(i-1)+1:2*i)*(FULLdata(j-s+2-i,1:2)'-mu);
                end						%end the if clause
            end							%end the iteration over the lagged y's.
        pi_for=yt_s(1,1);               %select the inflation forecast
        dgap_for=yt_s(2,1);
        piForecast_t_temp=[piForecast_t_temp;pi_for]; %stack them going down on forecast date
        dgForecast_t_temp=[dgForecast_t_temp;dgap_for]; %stack them going down on forecast date
    end
    piForecasts(1:N,t)=piForecast_t_temp;					%stack them horizontally
    dgForecasts(1:N,t)=dgForecast_t_temp;					%stack them horizontally
end
piForecasts=piForecasts(:,1:N-final);
dgForecasts=dgForecasts(:,1:N-final);
%display('1st column is forecasts of inflation in 2003:4 from different perspectives, and so on')
%display('first row is the actuals, while second row is one-period ahead, third row is two-periodsd ahead, etc.')
clear t j piForecast_t_temp dgForecast_t_temp s F_power constant yt_s i pi_for dgap_for mu k F_Forecasting F size_of_forecast_matrix

%STEP 4: Building one-step ahead forecasts for model
for j=1:upper
    power_lambda(j,1)=lambda*(1-lambda)^(j-1);
end
clear j
E=[];
for i=1:N-final             %iterate over calendar date strating from 2003:4 and moving backwards
    temp1=piForecasts(:,i)+alpha*dgForecasts(:,i);
    temp1=temp1(2:upper+1,1);
    temp2=dgForecasts(2,i)+gap(N-1,1);
    E=[E (alpha*lambda/(1-lambda))*temp2 + power_lambda'*temp1];
end
E=flipud(E');   %converting into a time series from 1996:1 to 2003:4
actual=FULLdata(final+1:N,1);
Comp=[E actual];
Comp = detrend(Comp,'constant');
modelMSE=1000*(sum((Comp(:,1)-Comp(:,2)).^2)/(N-final))
clear E

%STEP 5: Building one-step ahead forecasts for VAR
E2=piForecasts(2,:);
E2=flipud(E2');
Comp2=[E2 actual];
Comp2 = detrend(Comp2,'constant');
varMSE=1000*(sum((Comp2(:,1)-Comp2(:,2)).^2)/(N-final))
clear E2

%STEP 6: Building one-step ahead forecasts for AR
if final==183
    for i=1:N-final
        E3(1,i)=0.000872+0.536103*FULLdata(N-i,1)+0.371216*FULLdata(N-i-1,1);
    end
elseif final==167
    for i=1:N-final
        E3(1,i)=0.000995+0.531069*FULLdata(N-i,1)+0.369374*FULLdata(N-i-1,1);
    end
end
E3=flipud(E3');
Comp3=[E3 actual];
Comp3 = detrend(Comp3,'constant');
arMSE=1000*(sum((Comp3(:,1)-Comp3(:,2)).^2)/(N-final))

improvement=1-modelMSE/min([varMSE arMSE])